Self— forced gravitational waveforms for Extreme and Intermediate mass ratio inspirals 
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We present the first orbit-integrated self force effects on the gravitational waveform for an I(E)MRI 
source. We consider the quasi-circular motion of a particle in the spacetime of a Schwarzschild black 
hole and study the dependence of the dephasing of the corresponding gravitational waveforms due 
to ignoring the conservative piece of the self force. We calculate the cumulative dephasing of the 
waveforms and their overlap integral, and discuss the importance of the conservative piece of the self 
force in detection and parameter estimation. For long templates the inclusion of the conservative 
piece is crucial for gravitational-wave astronomy, yet may be ignored for short templates with little 
effect on detection rate. We then discuss the effect of the mass ratio and the start point of the 
motion on the dephasing. 
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I. INTRODUCTION AND SUMMARY 

The detection of gravitational waves (henceforth GW) 
and the onset of the new field of gravitational-wave as- 
tronomy is one of the most exciting challenges for sci- 
ence in the XXI century, completing what is sometimes 
alluded to as Einstein's Unfinished Symphony. The de- 
tection of GW will open a new window onto the uni- 
verse, that in addition to revealing exciting information 
on exotic systems such as black holes or cosmic strings is 
expected to also unravel as yet unexpected sources. 

One of the interesting sources for low-frequency GW 
are the so called I(E)MRI sources, or Intermediate (Ex- 
treme) Mass Ratio Inspirals. Those are the GW emitted 
by a system including a smaller compact object whose or- 
bit decays into a much larger massive black hole (MBH) . 
Typical sources are stellar mass black holes inspiraling 
into a supermassive black hole, like those residing at the 
center of galaxies, and also IMBHs (intermediate mass 
black holes) inspiraling into MBHs. The importance of 
such sources is that because of the extreme mass ratio 
the smaller compact object can be viewed as a test par- 
ticle, thus probing the spacetime of the larger black hole 
and its surroundings. Inter alia, such sources will al- 
low us to test directly the Kerr hypothesis, and allow us 
to map the spacetime surrounding such exotic objects. 
Moreover, the detection of I(E)MRIs will allow us to de- 
termine the mechanisms that shape stellar dynamics in 
galactic nuclei with unprecedented precision [l]. 

The orbits of I(E)MRIs are typically highly relativis- 
tic, and exhibit exciting phenomena, e.g. extreme perias- 
tron and orbital plane precessions. Because the orbital 
evolution time scale ("radiation reaction time scale") is 
much longer than the orbital period(s), over short time 
scales the orbit is approximately geodesic, yet on long 
time scales it deviates strongly from geodesic motion of 
the background. Instead, the smaller objects moves along 



a geodesic of a perturbed spacetime. Alternatively, one 
may construe the orbit as an accelerated, non-geodesic 
motion in the spacetime of the unperturbed central ob- 
ject, where the acceleration is caused by the self force 
(henceforth SF) of the smaller object @ ■ 

Detection and parameter estimation of GW from 
I(E)MRIs relies on the construction of theoretical tem- 
plates. A number of approximation schemes for such 
templates are available. Firstly, the energy balance ap- 
proach ("the radiative approximation") uses balance ar- 
guments for otherwise conserved quantities, and relates 
the flux in these quantities to infinity and down the event 
horizon of the black hole with the particle's orbit, so 
that the latter can be adjusted to agree with the fluxes 
H ■ As the orbital evolution time scale is typically much 
longer than the orbital period(s), the radiative approx- 
imation is very satisfactory during the adiabatic phase 
of the motion. As the particle's orbit is affected by the 
fluxes away from it, when the orbit is not stationary one 
encounters complex retardation effects. Most currently- 
available EMRI waveforms have been obtained by such an 
approach. This approach, however, ignores conservative 
effects that do not register in the constants of motion. 

These retardation effects are completely avoided when 
one considers a local approach to orbital evolution in 
terms of the SF. (One should bear in mind, however, that 
the SF itself is a non-local quantity, with contributions 
arising from the quasilocal neighborhood of the particle 
and possibly beyond [|[.) In addition, the local approach 
to the calculation of orbital evolution via the SF is not 
restricted to the adiabatic regime, it avoids the compli- 
cations associated with the rate of change of the Carter 
constant, and, most importantly, it includes also conser- 
vative effects that are discarded when one uses balance 
arguments. Over the last decade much progress has been 
made in the computation and understanding of the SF 
in General Relativity (for recent reviews of the self force 
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in General Relativity see @, Q). 

The computation of the fully relativistic SF allows 
one to include conservative effects in the waveform tem- 
plates, and study the importance of the conservative ef- 
fects. True self consistent orbit and waveforms include 
the instantaneous solution of the coupled SF integrated 
equations of motion and the perturbation equations, or 
equivalently the interaction of the particle with its own 
field over its half-infinite past world line @. Very re- 
cently, for the scalar field toy model, such self consistent 
Schwarzschild orbits and waveforms were presented [7|. 
Here, we are making the simplifying assumption that the 
effects of the difference between the SF that is calcu- 
lated for the actual orbit (the self consistent approach 
d, Q) and that which is calculated for a geodesic of the 
same instantaneous orbital parameters, is smaller than 
the effects of the latter and hence negligible at first or- 
der. This approximation is valid for as long as the orbital 
evolution is adiabatic, that is as long as the orbital evolu- 
tion time scale is much longer than the orbital period(s). 
In a Schwarzschild background of mass M, the adiabatic 
approximation holds when the mass ratio 77 := n/M is 
such that e ^> rj 1 / 2 , where e measures the distance to the 
innermost stable orbit, specifically e = p — 6 — 2e where 
pis the semilatus rectum and e is the orbital eccentricity 
[9(. In practice, our approximation is to a leading or- 
der in 77 beyond geodesic motion. We neglect terms that 
are linear in second-order SFs, although our method is 
amenable to their inclusion when they become available. 
This approximation is valid for at least a part of the rel- 
evant parameter space [10| . but as their inclusion would 
contribute linearly to the dephasing, the contribution of 
the conservative piece of the SF (hereafter CSF) may be 
isolated as is done here. Using true self consistent wave- 
forms will both produce more accurate waveforms, and 
allow us to test the accuracy of this approximation. Most 
importantly, our approach allows us to see for the first 
time the effect of the CSF on GW emitted from IMRI 
sources. 

We present here the first waveforms obtained with in- 
clusion of the CSF, and study its effect within the sim- 
ple class of quasi-circular orbits around a Schwarzschild 
black hole. Specifically, we study the effect of the sys- 
tem's mass ratio on the dephasing that occurs when one 
neglects the CSF. We find weak dependence of the de- 
phasing on the mass ratio, in accord with expectations 
based on the scaling of the number of orbits with the in- 
verse of the mass ratio, and the scaling of the dephasing 
effect of the CSF per orbit with the mass ratio. We also 
find that the dephasing depends quadratically on the ini- 
tial point of the motion for the range of parameters we 
tested. We reiterate that second-order dissipative effects 
are ignored in this Paper. Their inclusion will guarantee 
the full consistency of the model, and will be compa- 
rable to the self-consistent approach. The inclusion of 
the second-order dissipative effects awaits further devel- 
opment to both theory and computational techniques. 

The organization of this Paper is as follows: In Section 



HI! we discuss the computational and numerical methods 
that we use. In Section UTT1 we discuss our results for the 
orbits pil Al) . the waveforms ([III B[) and the dependence 
of the dephasing on the mass ratio and the initial point 
of the motion (|HIC[) . 



II. METHOD 

We use the fully relativistic SF obtained by Barack and 
Sago [ll| for circular Schwarzschild geodesies to drive the 
orbital evolution (our computation allows for an easy re- 
placement with a different force expression, say one that 
includes second order dissipative effects, or spin-orbit 
coupling effects when becoming available), and compare 
the resulting waveforms with those obtained from the 
energy balance approach and those obtained when only 
dissipation is left in the SF, setting by hand the CSF 
to vanish. In practice, we consider a point source \i in 
a quasi-circular orbit around a Schwarzschild black hole 
M with a mass ratio 77 :— (J./M, and the motion starts at 
a value of the semilatus rectum po until it gets close to 
the ISCO at p — 6. Such sources are relevant to the NGO 
capabilities: NGO will have the capability to detect GW 
emitted by an IMBH in the mass range 10 2 ~ 4 Mq spiral- 
ing into a MBH in the mass range 3 x 10 5 — 10 7 M© such 
that the mass ratio is 10~ 3 — I0~ 2 out to cosmological 
redshift z ~ 2 — 4. In addition, advanced LIGO could 
detect compact stellar sources spiraling into an interme- 
diate mass black hole (IMBH) in the same mass ratio 
range [12j. We specialize below to this mass ratio range, 
r) € [f0~ 3 — I0~ 2 ]. Although the linearized approach 
used here is intended to be used only when 77 <§; I and 
one may not simply extend the range of rj to high values 
and still expect accurate results, the error involved from 
neglecting 0(r) 2 ) terms in the self force (specifically its 
dissipative piece) is comparable to the accuracy of our 
computation. In this sense, to within the accuracy of our 
numerics, we are justified in studying the IMRI case as 
long as we do not raise r\ beyond 10 ~ 2 . 

We integrate the equations of motion using the 
Barack-Sago SF which was calculated for momentary cir- 
cular geodesies. As the Barack-Sago SF is tabulated for 
a select choice of orbital radii, we interpolate to inter- 
mediate values such that the original accuracy is main- 
tained. Specifically, we match two asymptotic expan- 
sions: at large distances (which we take in practice to be 
r > 8M) we take the standard post-Newtonian expan- 
sion for the luminosity in gravitational waves, and con- 
struct from it the temporal component of the SF. To 5.5 
PN order the PN expression does not provide us with suf- 
ficient accuracy to reproduce all the r > 8M data points 
of the Barack-Sago data. We therefore add an effective 
remainder term that appears like a 6PN term, and fit its 
two free parameters to agree with all the large distance 
tabulated data to all significant figures. The radial com- 
ponent, or CSF, is modeled by a PN-like expansion with 
four free parameters. At short distances we expand the 
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SF such that convergence is fast and only four free pa- 
rameters are needed for either component. Specifically, 
we expand the radial component about the ISCO, the In- 
nermost Stable Circular Orbit at p — 6. The expansion 
functions are as follows: 
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where PN5.5 stands for the standard -y— post-Newtonian 
expression (converting Eq. (3.1) in [13l | from luminosity 
to /' - see Appendix [S| . Fitting the free parameters, we 
find the values appearing in Table 1. 



TABLE I: The fit parameters for the self force. These param- 
eters reproduce the accuracy of [llj to all significant figures 
for all data points. Our results for a^ 6L are very inaccurate 
predictions for the corresponding PN parameters, as our fit 
ignores all higher-order terms. 
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The simplicity of our model allows us to easily sep- 
arate the CSF effects. Specifically, for quasi-circular 



Schwarzschild orbits we may write / 
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/i^ F + fr F $r wri ere the last term on the right hand 
side (RHS) is purely conservative, and the first two are 
purely dissipative. We may therefore study the conser- 
vative effects by turning off by hand the last term on the 
RHS. 

We integrate the SF driven orbit using two indepen- 
dent methods: one method is the osculating orbit ap- 
proach 15] [specifically Eqs. (43)-(47) therein], with spe- 
cial care given to the requirement that the orbit is quasi- 
circular. Specifically, free evolution may take the orbit 
away from quasi-circularity because integration using the 
osculating geodesies method cannot keep the value of the 
eccentricity as precisely zero. As both variables a and (3 
(see [lU for definitions) are dynamical, the eccentricity e 
must evolve along the orbit too. This behavior is shown 
in Fig. [T] Interestingly, the inclusion of the conservative 
piece of the self force amplifies the resulting eccentricity. 
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FIG. 1: The osculating-code variables a and j3 as functions 
of the time t (upper two panels), and the effective eccentricity 
e of the orbit (lower panel) as a function of t in the osculating 
case. We show the results from the osculating orbits method 
with the inclusion of the CSF (solid curves) and with turning 
off the CSF (dotted curves). In all cases shown here 77 = 10 -2 
and po = 10. 

The second method is the direct integration of the or- 
bit. The direct integration method takes the local equa- 
tion of motion to be u 13 V pu a = /g F ("Newton's sec- 
ond law," with covariant differentiation compatible with 
the background metric) and integrates its solution. Both 
codes are numerically stable and convergent. Specifically, 
the osculating code converges with 5 th order, and the di- 
rect code converges with 4 th order (Fig. [5]). 

Comparing the two independent methods for finding 
the SF driven waveforms is not trivial because of a diffi- 
culty in finding identical initial conditions. Specifically, 
the osculating method requires as initial data only the 
specification of the initial position vector, which in our 
case is taken to be a circular geodesic at some initial po- 
The direct method, however, requires both the position 
and the velocity vectors to be specified, such that the 
constraint equation u^u^ = —1 is also satisfied. The 
main difficulty is that in the initial data for the osculat- 



4 



% 




(B) 1 


1 







300 400 



200 300 




6 4 -° 



2000 3000 



4::iji) 
t/M 



5000 6000 



FIG. 2: Convergence tests for the two codes. Top panels: 
the 4-position (A) and 4-velocity (B) for the osculating code 
as functions of the azimuthal angle ip. Lower panels: the 
4-position (C) and 4-velocity (D) for the direct code and a 
function of time t. In all cases shown here r\ — 10~ 2 and 
r = 10M, and the CSF is included. 



ing method the initial u 7 = 0, such that it corresponds to 
an incorrect initial radial velocity for the direct method. 
In practice, we find the initial data for the direct code by 
generating the orbital parameters at po by running the 
osculating code from some p 3> po down to po, and then 
take the position and velocity vectors at po as the initial 
data for the direct method. The residual disagreement in 
the initial data can be controlled to be compatible with 
our numerical error tolerance. We then use the obtained 
orbits to generate the waveforms using a code for the 
sourced Teukolsky equation with hyperboloidal slicing, a 
code which converges at 2 nd order [14j . 

As noted above, our approximation holds only for as 
long as £ > pT 1 !" 1 . We therefore do not integrate the 
equations of motion in practice all the way down to the 
ISCO at p = 6, and stop the integration at a finite dis- 
tance from the ISCO. In practice, we stop the integration 
at Pfinai = 6.15 ± 0.10. Stopping at pgnai is enough to es- 
timate the dephasing at the ISCO: the dephasing A* is 
a smooth function of the time t along the particle's world 
line. In practice, we extrapolate r as a function of t to the 
ISCO to determine the time at which the particle arrives 
at the ISCO, and then we extrapolate the phase of the 
waveform to the same value of the time to estimate the 
phase of the waveform when the particle arrives at the 
ISCO. We may then find the difference of the total phases 
between two waveforms to find the dephasing A$. 



III. RESULTS 
A. The orbit 

We next choose rj — 10~ 2 and po = 10. There is of 
course nothing special about this choice of the parame- 



ters, except that we couldn't make a much higher choice 
for rj and justify it with the linearized approximation 
used. Below we study the dependence of the effect on 
the parameters rj, pq. The orbit is displayed in Figs. l3l4"l5l 
and [6] for the three codes. Notably, the two independent 
self force codes reproduce the orbit to high level of agree- 
ment, with a difference much smaller than the difference 
between either and the orbit generated in the energy bal- 
ance approach. This difference is attributed to the effect 
of the conservative piece of the self force. 
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FIG. 3: The Orbit. The 4-position for three orbital evolution 
codes: energy balance (dotted), direct evolution (dashed), 
and osculating code (solid). 
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FIG. 4: The Orbit. The 4-velocity for three orbital evolution 
codes: energy balance (dotted), direct evolution (dashed), 
and osculating code (solid). 

The orbit, of course, is a gauge dependent quantity. 
Indeed, the position vector changes trivially under gauge 
transformations, x^ — > x^ + £ M . We can, however, cre- 
ate gauge invariant quantities in a specific gauge choice 
(in our case, the Lorenz gauge), and then those quan- 
tities are guaranteed to remain unchanged in any other 
gauge. Two independent gauge invariant quantities are 
u t "g rav itational redshift" , "helical Killing vector of the 
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FIG. 5: The Orbit. The shape of the orbit for the three 
orbital evolution codes: energy balance (dotted), direct evo- 
lution (dashed), and osculating code (solid). 
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FIG. 6: The Orbit. The shape of the orbit for the three 
orbital evolution codes: energy balance (dotted), direct evo- 
lution (dashed), and osculating code (solid). A small portion 
of the orbit shown in Fig. [5] is magnified to show detail. 



perturbed spacetime" ) and the angular frequency [16j ■ 
In Fig. [7] we plot u* as a function of fl with and without 
the conservative piece of the self force. Notably, to the 
accuracy of our numerical computation the two curves 
overlap. That is, we find - as expected - that u l as a 
function of fi is insensitive to the conservative piece of 
the self force 0. This conclusion implies that when an 
actual data stream is used and this gauge invariant figure 
is plotted, one may use a simplified radiation-reaction 
scheme, that does not include the conservative effects in 
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FIG. 7: The gauge invariant figure of u as a function of 
the angular frequency fl. The curve without (dashed, o) and 
with the conservative effects (solid, □) are shown, together 
with equal- f increment marks starting at t = 1, 500 M (at the 
bottom left) in increments of 1, 000 M. The two curves are in- 
distinguishable to the numerical accuracy of our computation 
(using the osculating orbits code in both cases). 



its analysis. 

There is, however, an aspect of the gauge invariant fig- 
ure (J7J) that is sensitive to the conservative effects, specif- 
ically the speed with which the data point moves along 
the curve. The way the conservative effects are mani- 
fested in the gauge invariant plot is not is the shape of 
the curve, but in the time it takes the signal to move 
along it. One may therefore observe the CSF effect by 
monitoring the motion of the data point representing the 
system along its curve on the u* — O plane. We note that 
one additional effect of the CSF is the shift in the ISCO 
p^ . which we do not consider here. 



B. The waveforms 

We show the waveform for the case that the CSF is 
turned off in Fig. [5] for the same parameters discussed 
above, specifically po = 10 and rj = 10~ 2 . We find that 
the waveforms obtained with the SF keeping only its dis- 
sipative pieces (and turning off its CSF) overlap with 
the energy balance waveform. Notice that the two wave- 
forms in the figure are indistinguishable. The calculation 
method is very different in these two cases: In the SF case 
the orbital evolution is local; the orbit evolves because of 
a local force acting on the particle; in the energy balance 
case the fluxes to infinity and down the event horizon 
are calculated, and then the energy escaping over a pe- 
riod of the orbit is removed from the particle, and a new 
orbit with the new values of the constants of motion is 
found. The agreement of the waveforms is therefore a 
non-trivial test of the correctness of the calculation. The 
two waveforms overlap nearly exactly, with total cumu- 
lative dephasing at the order of 10 -3 radians. 
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FIG. 8: The real part of the waveform ip4 when the CSF 
is turned off (solid), and using the energy balance approach 
(dashed, o). The SF waveform was calculated with the orbit 
evolving using the osculating method [TH ]. 



We next reintroduce the CSF. We integrate the SF 
driven orbit using two independent methods: the oscu- 
lating orbit approach and direct integration of the orbit. 
Figure [9] shows the waveform for the energy balance ap- 
proach (same as in Fig. [3]) and for the two independent 
methods of calculating the SF driven orbit (including 
the CSF). The latter two waveforms are in agreement 
with each other with small dephasing (see below) be- 
tween them, and a much larger dephasing of either with 
the energy balance waveform. 
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FIG. 9: The real part of the waveform tpi when the CSF is 
included (solid for the osculating orbit method, and dashed for 
the direct integration method), and using the energy balance 
approach (dotted). The three waveforms are in phase at the 
beginning of the waveforms, but the former two dephase with 
time from the latter (and to a much smaller extent from each 
other). 

The waveform dephasing of Fig. [5] is shown explicitly 
in Fig. [TO] We find that the total cumulative dephasing 
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FIG. 10: The dephasing of the waveforms for the same data 
shown in Fig. [9] We show the dephasing from the energy bal- 
ance waveform of the waveforms obtained from both methods 
to calculate the orbital evolution when the conservative effects 
are included, and the dephasing of the two latter methods 
from each other. Dashed curve: dephasing of the waveforms 
from the osculating orbit method from the energy balance 
method; Dotted curve: dephasing of the waveforms from the 
direct integration method from the energy balance method; 
Thick solid curve: the average dephasing of the waveform. 
Dash-dotted curve: The difference between the dotted and 
dashed curves. Insert: same as in the main figure, with the 
extrapolated dephasing down to the ISCO shown in a thin 
solid curve. 



of the waveforms at the endpoint of the evolution at pfinai 
is A</> = 10.3 ± 0.4 radians. The error estimate comes 
from the numerical errors in each calculation method and 
from residual incompatibility of initial data in the two 
SF driven cases, which we estimate by comparing the 
waveforms obtained from the osculating orbit approach 
and the direct integration approach. 

The dependence of A<E> on the position p is a simple 
monotonic unction of the time (Fig. I10|) which we can 
extrapolate from pgnal down to the ISCO at p = 6. We 
find that at the ISCO the dephasing is A$ = 14 ± 1 
radians. Dephasing of 14±1 radians corresponds to about 
2.2 cycles over the entire motion of the particle over 107.8 
cycles. We next consider the following simulation of a 
detection event. Say the actual data stream is modeled 
by the waveforms obtained with the full SF expression, 
i.e., including the CSF, and that the theoretical template 
is obtained by turning off the CSF, or equivalently by 
using the energy balance waveform. By how much is the 
overlap integral of the waveforms reduced because of our 
ignorance of the CSF? 

Figure [TT] shows the overlap integral as a function of 
the time, when we take a window of length L of the en- 
ergy balance waveforms (specifically from their late chirp 
part), and shift it along the waveform of the full SF. At 
each point we calculate the overlap integral of the win- 
dow with a local piece of the second waveform, and plot 
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FIG. 11: The local overlap integral for a window of length 
L = 350M (dash-dotted) and of length L = 820Af (solid) 
taken from the energy balance waveform and a piece of equal 
length of the full SF (including conservative effects) waveform 
obtained from the osculating orbit method, as a function of 
the starting point of the latter. The window is taken here 
from the end of the waveform (the late chirp part). The global 
maximum of the overlap integral is 0.9900 (for L — 350M) 
and 0.9594 (for L = 820M). 



FIG. 12: The maximal overlap integral for windows of varying 
length L taken from the energy balance waveform and a piece 
of equal length of the full SF (including conservative effects) 
obtained from the osculating orbit method, as a function of 
the starting point of the latter. The circles are data points, 
and the curve is the quadratic fit C max = 1.0010 — 1.7988 x 
10" 5 L - 3.9154 x 10" 8 L 2 , with fit parameter R 2 = 0.9983. 



C. Varying the values of parameters 



the local overlap integral as a function of the start time 
of the window. Because the window was taken from the 
late part of the waveform, we find that the overlap in- 
tegral is very small at first, and becomes large only at 
the late part of the other waveform. We then take the 
maximum of the overlap integral to be the one corre- 
sponding to the chosen window L. More precisely, we 

calculate C max = max T <^ (*W ( t ~ T )> _ 

^<V> SF (t)IV> SF «><V> EB (0IV' EB (*)> 
Larger values of L would reduce the overlap integral even 
further. In Fig. [T^] we show C max as a function of the win- 
dow size L. As L increases, at some value C max would 
drop below a pre-determined value that marks our toler- 
ance for detection or parameter estimation. Many times 
this threshold is taken to be C = 0.96, because then 
detection rate would drop by 10%. Here, this thresh- 
old is obtained when ^threshold = 816.6 M, which corre- 
sponds to just over 14 wavelengths of the emitted GW. 
If L > ^threshold) the exclusion of the CSF effects would 
cause a significant drop in the C max that would reduce the 
detection rate by 10% or more. In such a case ignoring 
the CSF would have an important effect on detection or 
parameter estimation of the GW. However, short wave- 
forms (i.e., L < -^threshold) do not require the CSF effects 
to be included if reduction of the detection rate by less 
than this tolerance is acceptable. The overlap integral 
increases rapidly as the template window is taken from 
earlier times. This result suggests that when other parts 
than the very end of the waveform is of interest, the full 
SF is even less significant than we have found. 



Our model of quasi-circular Schwarzschild orbits de- 
pends on two variables: the mass ratio 77 and the start 
point pq. Here we vary each parameter independently 
and find the dependence of the dephasing A$ [between 
the osculating orbits case (that includes the CSF) and the 
energy balance case (that neglects the CSF and considers 
only dissipative effects)] on either parameter. 



1. Varying the mass ratio rj 

First we study the variation of the dephasing with 
changing the parameter 77, the mass ratio. The greatest 
problem with varying 77 is its effect on the computation 
time. On the one hand we cannot justifiably increase rj 
beyond 10 -2 , because then the linearization approxima- 
tion breaks down. On the other hand, lowering 77 to very 
small values, while satisfying the linearization require- 
ment more confidently, results in longer physical evolu- 
tion times and correspondingly also longer computation 
times. 

In practice we reduce the value of 77 by a full order of 
magnitude, and sample values in the range [10~ 3 — 10~ 2 ], 
and fix po — 8. (We decrease the value of po from its 
previous value of 10 to save on computation time for the 
lower values of 77.) 

Figure [T3] shows a family of curves displaying the 
dephasing between the cases of the osculating orbits 
method and the energy balance method. Each curve in 
Fig. [13] ends at the point we stop the integration when 
the orbits gets too close to the ISCO for the adiabatic ap- 
proximation to still hold. These curves are very smooth 
and simple functions, which we extrapolate to the times 
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FIG. 13: The dephasing A<& as function of time for a family 
of mass ratios r\ for the values (from right to left): r\ =0.001, 
0.002, 0.003, 0.005, 0.006, 0.007, and 0.010 (solid curves). 
The circles display the extrapolated values when the particle 
arrives at the ISCO, and the dashed line shows their average 
at 1.726 ± 0.021. In all cases the motion starts at po = 8. 



at which the particle arrives at the ISCO. The extrapo- 
lated value of the dephasing at the ISCO is also shown 
in Fig. Q2J We find only little variation in the dephasing 
as r\ changes over a full order of magnitude, consistent 
with the expectation that the dephasing has only a weak 
dependence on r\. Indeed, one could expect from scal- 
ing arguments that the dephasing per orbit scales with r\ 
while the number of orbits scales with r\~ x (both scalings 
are indeed found in our simulations — see Fig. I14[) . so 
that the total dephasing is at the leading order at 0(rp). 
Here we show that not only is the dephasing at 0(1), but 
in the range tested and with our numerical resolution is 
indistinguishable from a constant value, or at the most 
is a very weak function of r\. 

In Fig. Q3] we show five dimensionless quantities con- 
structed from the evolution time T from po down to the 
ISCO, the total phase of the waveform in the energy bal- 
ance and osculating methods case (including the CSF 
term) , and the differences in arrival time and the dephas- 
ing between the latter two. We find that T ~ 0(?7 _1 ), 
that A$ - 0{i] ), AT/M ~ 0{rf) and that in both 
the energy balance and osculating orbits cases the total 
phase $ is a linear function of rj with a very small slope, 
that is at the magnitude of our computational error. We 
therefore cannot rule out that we see in addition to the 
leading 0(?7 _1 ) term also a higher-order O(rj ) term. No- 
tice that the last four quantities are comparable to each 
other, and that the variation in all five (at the most 10%) 
is very small compared with the full order of magnitude 
variation in rj. 



2. Varying the initial position po 

Next we fix the mass ratio r\ and vary the starting 
point of the motion pq. In practice we choose the value 



FIG. 14: Five dimensionless quantities (denotes collectively 
by iV) as functions of the mass ratio n: fiT/M 2 (divided by 10 
to keep the scale similar with the other four quantities) (o), 
the dephasing A4> (*), the difference in arrival time AT/M, 
and »7<I> for the energy balance case (□) and the osculating 
orbits case (o). Here, A<£> and AT/M are between the oscu- 
lating orbits and the energy balance cases, and T is the total 
time of motion from po down to the ISCO. In all cases the mo- 
tion starts at po = 8 and the values shown are extrapolations 
to the ISCO. 



a. s 




FIG. 15: The dephasing A3> as function of the initial posi- 
tion po. The circles display the extrapolated values when the 
particle arrives at the ISCO, and the curve shows a quadratic 
fit. In all cases the mass ratio n = 10~ 2 . 



of rj = 10~ 2 . Figure fl5l shows the dephasing between the 
osculating orbits case (that includes the CSF) and the 
energy balance case (that neglects the CSF) for the range 
Po € [8,10]. Naturally, the dephasing grows with po- The 
data presented are consistent a quadratic dependence of 
A<E> on po- 

After the completion of this work we became aware of 
Warburton et al [Hj]. The approach of 18 1 is similar to 
our osculating method, except that [l8| estimates the de- 
phasing by the difference in the azimuthal angle ip of the 
orbit, whereas we actually compute the waveforms and 
find their dephasing. The generalization of our quasi- 
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circular orbit to bound orbits of varying eccentricity — 
as is done in [l8[ — is straightforward, as the osculat- 
ing orbits equations of motion already include the ec- 
centricity parameter. Lastly, we compute the waveforms 
using two independent computational methods, specifi- 
cally the osculating method and the direct method. The 
direct method does not appear to us to be convenient for 
generalization to generic Kerr orbits, and the osculating 
method has a clear advantage over it for such orbits. 
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Appendix A: The 5.5 Post— Newtonian term 

The PN 5 . 5 term used in Eq. © is given by (43]): 
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